MODULE smm
  USE GLOBVAR
  USE NRUTIL
  USE PROBABILITY
  USE MATRIX
  USE statistics
  USE random
  IMPLICIT NONE
  
    CONTAINS
   REAL(8) FUNCTION mom(thetain) 
   !SUBROUTINE mom(theta_ini,beta_f,check,N_int,sc)
        IMPLICIT NONE
        REAL(8), INTENT(IN):: thetain(:) !initial guess

        REAL(8):: Ls(N_inds,N_testf),Ep(N_inds,N_testf),as,pen,pen2
        REAL(8):: L0s(N_inds),eta(N_inds,N_testf),e_L0(n_inds), e_eta(n_inds),mom_simu(N_mom),dif(N_mom,1)!, latent(N_inds*N_testl,N_test),ps(N_testl,N_test)
        REAL(8):: inters(n_inds,n_levelf,6,3),interfs(n_inds,n_levelf,3)
        REAL(8),ALLOCATABLE:: AVSE(:),cors(:,:)
        INTEGER:: dufs(n_inds,n_levelf),dus(n_levelf),inds(N_inds,n_levelf) 
        ! parameters
        REAL(8):: beta_mu(6,1),beta_eta(n_x_eta,1),var_l,fbar(N_testf),var_e(N_testf), etav,delta(N_levelf),beta_level(n_l,1),beta_y(1:2,1)
        INTEGER:: Ds(N_inds,N_testf),indexf(n_testf,2), artfs(n_inds,n_levelf,6) 
        INTEGER:: i,j,k,jj,kk,m,t,n,l,tt,mmm
        
        !!!!!!!!!!!!!!!!!!Simulate Moments
        CALL SET_SEED(1438914432,727629,19890414,111111)
        beta_mu(1:n_x_l0,1)=thetain(1:n_x_l0)                                              ! Initial Latent skill
        DO i=2,4
          beta_mu(i,1)=1.40d0**(thetain(i))/(10.0d0*(1.0d0+1.40d0**(thetain(i))))
        END DO
        beta_mu(5,1)=1.0d0-1.40d0**(thetain(5))
        
        var_l=4.0d0*(1.4d0**thetain(n_x_l0+1)/(1.0d0+1.4d0**thetain(n_x_l0+1)))     ! variance of the shoch of initial latent skill
                fbar=1.0d0
        fbar(2:N_levelf) =thetain(n_x_l0+2:n_x_l0+N_levelf)                ! mimimum requriement for each difficulty level
        var_e=1.0d0
        var_e(2:N_levelf)=1.4d0**thetain(1+n_x_l0+N_levelf:n_x_l0+2*N_levelf-1)         ! variance of the shock at each difficulty level
        beta_eta(1:n_x_eta,1)=thetain(n_x_l0+2*N_levelf:n_x_l0+2*N_levelf+n_x_eta-1)           ! x for eta
        etav=1.40d0**thetain(n_x_l0+2*N_levelf+n_x_eta)                           ! variance of shock for eta
        delta=1.0d0
        DO i=3,N_levelf
                        IF (i==3) THEN
                              delta(i)=4.0D0 +thetain(n_x_l0+2*N_levelf+n_x_eta+i-2)
                         ELSE IF (i==6) THEN
                            delta(i)=4.0D0 +thetain(n_x_l0+2*N_levelf+n_x_eta+i-2)
                         ELSE
                             delta(i)=4*(1.40d0**thetain(n_x_l0+2*N_levelf+n_x_eta+i-2))/(1.0d0+1.40d0**thetain(n_x_l0+2*N_levelf+n_x_eta+i-2))  ! delta for each level
                         END IF
        END DO
        beta_level=0.0d0
        beta_level(n_levelf+1:n_l,1)=thetain(n_x_l0+2*N_levelf+n_x_eta+n_levelf-2+1:n_x_l0+2*N_levelf+n_x_eta+n_levelf-2+n_levelf)
        beta_y=0.0d0
        beta_y(1:2,1)=thetain(n_x_l0+2*N_levelf+n_x_eta+n_levelf-2+n_levelf+1:n_x_l0+2*N_levelf+n_x_eta+n_levelf-2+n_levelf+2)



        ! generate shock terms
        
        DO i=1,N_inds
            e_L0(i) =Sample_NORMAL(0.0d0,var_l)
            k=i-INT((i-1)/N_IND)*N_ind
            e_eta(i)=Sample_Normal(intermean(k),etav)
        END DO
        
        ! GENERATE INITIAL LATENT VALUES AND 
        DO i=1,N_inds
            k=i-INT((i-1)/N_IND)*N_ind
            L0s(i)=1.40d0**(DOT_PRODUCT(L0(k,:),beta_mu(:,1))+e_L0(i))
        END DO
        
        
       ! SIMULATE LATENT SKILL PROCESS ! HERE ONLY FOR LANGUAGE SKILLS 
        
        DO i=1,N_inds
            k=i-INT((i-1)/N_IND)*N_ind
            DO j=1,N_testf
                IF (j<firstf(k,1)) THEN
                    Ls(i,j)=-99.0d0
                ELSE IF (j==firstf(k,1) .and. m_f(k,j)==1) THEN
                    IF (x_testf(k,j,1)<-90) THEN
                        eta(i,j)=1.40d0**(e_eta(i))/((1.0d0+1.40d0**(e_eta(i))))
                    ELSE
                        eta(i,j)=1.40d0**(DOT_PRODUCT(x_testf(k,j,:),beta_eta(:,1))+e_eta(i))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testf(k,j,:),beta_eta(:,1))+e_eta(i))))
                    END IF
                    
                    kk=levelf(j)
68                  Ls(i,j)=L0s(i)*(exp(delta(kk)*eta(i,j)+Sample_Normal(0.0d0,var_e(kk))+beta_y(1,1)*monthf(j)+beta_y(2,1)*DBLE(timef(j))))
                    !IF (Ls(i,j)<0.0d0) GO TO 68 
                ELSE IF (j==firstf(k,1) .and. m_f(k,j)==-1) THEN
                    Ls(i,j)=-99.0d0
                ELSE IF (j>firstf(k,1) .and. Ls(i,j-1)<0.0d0 .and. m_f(k,j)==-1) THEN
                    Ls(i,j)=-99.0d0
                ELSE IF (j>firstf(k,1) .and. Ls(i,j-1)<0.0d0 .and. m_f(k,j)==1) THEN
                    IF (x_testf(k,j,1)<-90) THEN
                        eta(i,j)=1.40d0**(e_eta(i))/((1.0d0+1.40d0**(e_eta(i))))
                    ELSE
                        eta(i,j)=1.40d0**(DOT_PRODUCT(x_testf(k,j,:),beta_eta(:,1))+e_eta(i))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testf(k,j,:),beta_eta(:,1))+e_eta(i))))
                    END IF
                    kk=levelf(j)
77                  Ls(i,j)=L0s(i)*(exp(delta(kk)*eta(i,j)+Sample_Normal(0.0d0,var_e(kk))+beta_y(1,1)*monthf(j)+beta_y(2,1)*DBLE(timef(j))))
                    !IF (Ls(i,j)<0.0d0) GO TO 77
                ELSE IF (j>firstf(k,1) .and. Ls(i,j-1)>=0.0d0 .and. m_f(k,j)==-1 .and. (levelf(j)==levelf(j-1))) THEN
                    Ls(i,j)=Ls(i,j-1)
                ELSE IF (j>firstf(k,1) .and. Ls(i,j-1)>=0.0d0 .and. m_f(k,j)==-1 .and. (levelf(j) .ne. levelf(j-1))) THEN
                    kk=levelf(j)
                    IF (kk>1) THEN
                       Ls(i,j)=(beta_level(kk,1)*1.0d0+beta_level(n_levelf+kk,1)*Ls(i,j-1))  
                    ELSE 
                        Ls(i,j)=Ls(i,j-1)
                    END IF
                ELSE IF (j>firstf(k,1) .and. Ls(i,j-1)>=0.0d0 .and. m_f(k,j)==1 .and. (levelf(j)==levelf(j-1))) THEN
                    IF (x_testf(k,j,1)<-90) THEN
                        eta(i,j)=1.40d0**(e_eta(i))/((1.0d0+1.40d0**(e_eta(i))))
                    ELSE
                        eta(i,j)=1.40d0**(DOT_PRODUCT(x_testf(k,j,:),beta_eta(:,1))+e_eta(i))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testf(k,j,:),beta_eta(:,1))+e_eta(i))))
                    END IF
                    kk=levelf(j)
84                  Ls(i,j)=Ls(i,j-1)*(exp(delta(kk)*eta(i,j)+Sample_Normal(0.0d0,var_e(kk))+beta_y(1,1)*monthf(j)+beta_y(2,1)*DBLE(timef(j))))
                    !IF (Ls(i,j)<0.0d0) GO TO 84
                ELSE IF (j>firstf(k,1) .and. Ls(i,j-1)>=0.0d0 .and. m_f(k,j)==1 .and. (levelf(j) .ne. levelf(j-1))) THEN
                    IF (x_testf(k,j,1)<-90) THEN
                        eta(i,j)=1.40d0**(e_eta(i))/((1.0d0+1.40d0**(e_eta(i))))
                    ELSE
                        eta(i,j)=1.40d0**(DOT_PRODUCT(x_testf(k,j,:),beta_eta(:,1))+e_eta(i))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testf(k,j,:),beta_eta(:,1))+e_eta(i))))
                    END IF
                    kk=levelf(j)
92                  Ls(i,j)=((beta_level(kk,1)*1.0d0+beta_level(n_levelf+kk,1)*Ls(i,j-1)))*(exp(delta(kk)*eta(i,j)+Sample_Normal(0.0d0,var_e(kk))+beta_y(1,1)*monthf(j)+beta_y(2,1)*DBLE(timef(j))))
                    !IF (Ls(i,j)<0.0d0) GO TO 92    
                END IF
            END DO
        END DO
        
        ! GENERATE TASK PERFORMANCE !
        
        DO i=1,N_inds
            k=i-INT((i-1)/N_IND)*N_ind
            DO j=1,N_testf
                kk=levelf(j)
                IF (m_f(k,j)==-1) THEN
                    DS(i,j)=-99
                ELSE IF (Ls(i,j)>=SUM(fbar(1:kk)) .and. m_f(k,j)==1) THEN
                    DS(i,j)=1
                ELSE IF (Ls(i,j)<SUM(fbar(1:kk)) .and. m_f(k,j)==1) THEN
                    DS(i,j)=0
                END IF
            END DO
        END DO
                
         !!!!!!!!!!!!!!!!!From here !!!!!!!!!!!!!!!!!!!!!!!
        
    !!! cognitive tasks!!!!
    indexf=0
    DO i=1,N_inds
      DO j=1,N_testf
         IF (ds(i,j)>-90) THEN
             indexf(j,1)=indexf(j,1)+1
             IF (ds(i,j)>0.5) indexf(j,2)=indexf(j,2)+1
         END IF         
      END DO
    END DO
    
    DO i=1,N_testf
        IF (indexf(i,1)>0.5) THEN
            mom_simu(i)=DBLE(indexf(i,2))/DBLE(indexf(i,1)) ! PASSING RATE FOR EACH Cognitive TASK
        ELSE 
            mom_simu(i)=-99
        END IF
    END DO
    
    !!!! generate artc variables !!!!
    artfs=-1
    DO i=1,N_inds
        DO j=1,N_testf
            IF (ds(i,j)>-90) THEN
                kk=levelf(j)
                k=i-INT((i-1)/N_IND)*N_ind
                m=repf(k,j)
                artfs(i,kk,m)=ds(i,j)
                inters(i,kk,m,1:3)=x_testf(k,j,1:3)
            END IF
        END DO
    END DO
    
    !!!! passing rate at each level!!!!

    DO j=N_testf+1,N_testf+N_levelf
       m=0
       n=0 
       DO i=1,N_inds
           DO k=1,6
               IF (artfs(i,j-N_testf,k)>-0.5) THEN
                   m=m+1
                   IF (artfs(i,j-N_testf,k)==1) n=n+1
               END IF
           END DO 
       END DO
       IF (m>0) THEN
           mom_simu(j)=DBLE(n)/DBLE(m)
       ELSE
           mom_simu(j)=-99
       END IF
    END DO
    
    !!!! correlation across levels !!!!
    
    mmm=N_testf+N_levelf
    
   l=0 
   DO j=1,N_levelf-1
       IF (j==1) THEN
           tt=5
           kk=2
       ELSE IF (j==2) THEN
           tt=2
           kk=5
       ELSE IF (j==3) THEN
           tt=5
           kk=5
       ELSE IF (j==4) THEN
           tt=5
           kk=3
       ELSE IF(j==5) THEN
           tt=3
           kk=4
       ELSE IF (j==6) THEN
           tt=4
           kk=3
       END IF
      DO t=1,tt
        DO k=1,kk
          m=0
          n=0
          l=l+1
          DO i=1,N_inds 
             IF (artfs(i,j,t)==1 .and. artfs(i,j+1,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artfs(i,j,t)==1 .and. artfs(i,j+1,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(mmm+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(mmm+l)=-99
          END IF
        END DO
      END DO
   END DO
    
    mmm=N_testf+N_levelf+84
    
    l=0
    DO j=1,N_levelf
        
       IF (j==1) THEN
           tt=4
           kk=5
       ELSE IF (j==2) THEN
           tt=1
           kk=2
       ELSE IF (j==3) THEN
           tt=4
           kk=5
       ELSE IF (j==4) THEN
           tt=4
           kk=5
       ELSE IF(j==5) THEN
           tt=2
           kk=3
       ELSE IF (j==6) THEN
           tt=3
           kk=4
       ELSE IF (j==7) THEN
           tt=2
           kk=3
       END IF
        
      DO t=1,tt
        DO k=t+1,kk
          m=0
          n=0
          l=l+1
          DO i=1,N_inds 
             IF (artfs(i,j,t)==1 .and. artfs(i,j,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artfs(i,j,t)==1 .and. artfs(i,j,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(mmm+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(mmm+l)=-99
          END IF
        END DO
      END DO
    END DO
    
    mmm=N_testf+N_levelf+84+43
    l=0
    DO j=1,N_levelf-2
        IF (j==1) THEN
            tt=5
            kk=5
        ELSE IF (j==2) THEN
            tt=2
            kk=5
        ELSE IF (j==3) THEN
            tt=5
            kk=3
        ELSE IF (j==4) THEN
            tt=5
            kk=4
        ELSE IF (j==5) THEN
            tt=3
            kk=3
        END IF
        
      DO t=1,tt
        DO k=1,kk
          m=0
          n=0 
          l=l+1
          DO i=1,N_inds 
             IF (artfs(i,j,t)==1 .and. artfs(i,j+2,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artfs(i,j,t)==1 .and. artfs(i,j+2,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(mmm+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(mmm+l)=-99
          END IF
        END DO
      END DO
    END DO  
    
     mmm=N_testf+N_levelf+84+43+79
    
    
    indexf=0
    DO i=1,N_inds
         k=i-INT((i-1)/N_IND)*N_ind
      DO j=1,N_testf
          as=firstf(k,2)-monthf(j)-DBLE(timef(j))/4.0d0
         IF (ds(i,j)>-90 .and. as>-1.25d0) THEN
             indexf(j,1)=indexf(j,1)+1
             IF (ds(i,j)>0.5) indexf(j,2)=indexf(j,2)+1
         END IF
      END DO
    END DO
    
    DO i=1,N_testf
        IF (indexf(i,1)>0.5) THEN
            mom_simu(mmm+i)=DBLE(indexf(i,2))/DBLE(indexf(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
        ELSE 
            mom_simu(mmm+i)=-99
        END IF
    END DO
    
    mmm=2*N_testf+N_levelf+84+43+79
    
    indexf=0
    DO i=1,N_inds
         k=i-INT((i-1)/N_IND)*N_ind
      DO j=1,N_testf
          as=DBLE(firstf(k,2)-monthf(j)-DBLE(timef(j))/4.0d0)
         IF (ds(i,j)>-90 .and. as<-1.25d0) THEN
             indexf(j,1)=indexf(j,1)+1
             IF (ds(i,j)>0.5) indexf(j,2)=indexf(j,2)+1
         END IF
      END DO
    END DO

    
    
    DO i=1,N_testf
        IF (indexf(i,1)>0.5) THEN
            mom_simu(mmm+i)=DBLE(indexf(i,2))/DBLE(indexf(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
        ELSE 
            mom_simu(mmm+i)=-99
        END IF
    END DO
    
    mmm=3*N_testf+N_levelf+84+43+79
    
    l=0
    DO j=1,N_levelf-1
      DO t=1,5
          m=0
          n=0
          l=l+1
          DO i=1,N_inds 
             IF (artfs(i,j,t)>-0.5) THEN   
                 m=m+1
                 IF (artfs(i,j,t)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(mmm+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(mmm+l)=-99
          END IF
      END DO
    END DO
    mmm=3*N_testf+N_levelf+84+43+79+5*(N_levelf-1)  
    
    DUfs=-99
    inds=-99
    DO i=1,N_inds
        DO k=1,N_levelf
            DO j=1,6
                IF (artfs(i,k,j)==1 .and. j==1) THEN
                    dufs(i,k)=1
                    inds(i,k)=1
                ELSE IF (artfs(i,k,j)==0 .and. j==1) THEN
                    dufs(i,k)=j
                    inds(i,k)=0
                ELSE IF (j>1 ) THEN
                    IF( artfs(i,k,j-1)==0 .and. artfs(i,k,j)==0 .and. inds(i,k)==0) THEN
                       dufs(i,k)=j
                       inds(i,k)=0
                       IF (j<6 ) THEN
                           IF (artfs(i,k,j+1)==-1) THEN
                           dufs(i,k)=j+1
                           inds(i,k)=1
                           END IF
                       END IF  
                    ELSE IF (j>1 .and. artfs(i,k,j-1)==0 .and. artfs(i,k,j)==0 .and. inds(i,k)==1) THEN
                        dufs(i,k)=dufs(i,k)
                        inds(i,k)=1
                    ELSE IF (j>1 .and. artfs(i,k,j-1)==0 .and. artfs(i,k,j)==1 .and. inds(i,k)==0) THEN
                        dufs(i,k)=j
                        inds(i,k)=1
                    ELSE IF (j>1 .and. artfs(i,k,j-1)==0 .and. artfs(i,k,j)==1 .and. inds(i,k)==1) THEN
                        dufs(i,k)=dufs(i,k)
                        inds(i,k)=1   
                    END IF
                END IF
                
            END DO           
                
        END DO
    END DO
    
    DO k=1,n_levelf
        m=0
        n=0
        DO i=1,n_inds
            IF (dufs(i,k)>-90) THEN
               m=m+1
               n=n+dufs(i,k)
            END IF
        END DO
        IF (m>0) THEN
          mom_simu(mmm+k)=DBLE(n)/DBLE(m)
        ELSE 
           mom_simu(mmm+k)=-99
        END IF
    END DO
    
    mmm=3*N_testf+N_levelf+84+43+79+5*(N_levelf-1)+n_levelf  
    
        dus=0
    DO k=1,n_levelf
        DO i=1,N_inds
            IF (dufs(i,k)>-90) THEN
                dus(k)=dus(k)+1
            END IF
        END DO
    END DO
    DO k=1,n_levelf
        allocate(avse(dus(k)))
        m=0
        DO i=1,n_inds
            IF (dufs(i,k)>-90) THEN
                m=m+1
                avse(m)=dufs(i,k)
            END IF
        END DO
        mom_simu(mmm+k)=STDEV_V(AVSE)
        DEALLOCATE(avse)
    END DO
     mmm=3*N_testf+N_levelf+84+43+79+5*(N_levelf-1)+2*n_levelf 

    INDs=0
    interfs=0.0d0
    DO k=1,N_levelf
        DO i=1,N_inds
          DO j=1,6
              IF (inters(i,k,j,1)>-90) THEN
                  inds(i,k)=inds(i,k)+1
                  interfs(i,k,1)=interfs(i,k,1)+inters(i,k,j,1)
                  interfs(i,k,2)=interfs(i,k,2)+inters(i,k,j,2)
                  interfs(i,k,3)=interfs(i,k,3)+inters(i,k,j,3)
              END IF
          END DO
          IF (inds(i,k)>0) THEN
               interfs(i,k,1)=DBLE(interfs(i,k,1))/DBLE(inds(i,k))
               interfs(i,k,2)=DBLE(interfs(i,k,2))/DBLE(inds(i,k))
               interfs(i,k,3)=DBLE(interfs(i,k,3))/DBLE(inds(i,k))
          END IF
        END DO
    END DO
    
    dus=0
    DO k=1,n_levelf
        DO i=1,N_inds
            IF (inds(i,k)>0 .AND. dufs(i,k)>-90) THEN
            dus(k)=dus(k)+1          
            END IF
        END DO
    END DO
    
    DO k=1,n_levelf
        allocate(cors(dus(k),4))
        m=0
        DO i=1,n_inds
             IF (inds(i,k)>0 .AND. dufs(i,k)>-90) THEN
                m=m+1
                cors(m,1)=interfs(i,k,1)
                cors(m,2)=interfs(i,k,2)
                cors(m,3)=interfs(i,k,3)
                cors(m,4)=dufs(i,k)
             END IF
        END DO
        !IF (k>1) THEN
        mom_simu(mmm+3*(k-1)+1)=Correlation(cors(:,1),cors(:,4))
        mom_simu(mmm+3*(k-1)+2)=Correlation(cors(:,2),cors(:,4))
        mom_simu(mmm+3*(k-1)+3)=Correlation(cors(:,3),cors(:,4))
        !ELSE 
        !    mom_simu(mmm+3*(k-1)+1)=-99
        !    mom_simu(mmm+3*(k-1)+2)=-99
        !    mom_simu(mmm+3*(k-1)+3)=-99
        !END IF
        
            
        DEALLOCATE(cors)
    END DO
    

    
    
    !    
    !     !!!!!!!!!!!!! Now finish simulate the fake data and start to generat the moments !!!!!!!!!!!!!!!!!!!!!!!!!
    !     
    
    DO i=1,N_mom
        IF (ABS(moment(i))<1.0d-10) THEN
            dif(i,1)=((mom_simu(i)-moment(i)))**2.0D0
        ELSE 
            dif(i,1)=((DBLE((mom_simu(i)-moment(i))))**2.0D0)/DBLE(ABS(moment(i)))
        END IF 

    END DO
    
    mom=SUM(dif)
    
    
    
                OPEN(6424,file="dif.out")
                  DO i = 1,N_mom
                    WRITE(6424,'(200F16.8)') SQRT(dif(i,1)),mom_simu(i),moment(i)
                  END DO
                CLOSE(6424)
            
   
   END FUNCTION mom
    
       
    
END MODULE smm 